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For a previously published study of the titanium hep (a) to omega (lu) transformation, a tight- 
binding model was developed for titanium that accurately reproduces the structural energies and 
electron eigenvalues from all-electron density-functional calculations. We use a fitting method that 
matches the correctly symmetrized wavefuctions of the tight-binding model to those of the density- 
functional calculations at high symmetry points. The structural energies, elastic constants, phonon 
spectra, and point-defect energies predicted by our tight-binding model agree with density-functional 
calculations and experiment. In addition, a modification to the functional form is implemented to 
overcome the "collapse problem" of tight-binding, necessary for phase transformation studies and 
molecular dynamics simulations. The accuracy, transferability and efficiency of the model makes it 
particularly well suited to understanding structural transformations in titanium. 

PACS numbers: 71.15.Nc, 61.72.Ji, 63.20.-e, 62.20.Dc 



I. INTRODUCTION 

Titanium is a useful starting material for many struc- 
tural alloysjA however, the formation of the high-pressure 
omega phase is known to lower toughness and ductility. 2 
The atomistic mechanism of the transformation from the 
room temperature a phase (hep) to the high-pressure ui 
was recently elucidated by Ref. • The explication of the 
a— >lu atomistic transformation relied on the comparison 
of approximate energy barriers for nearly 1000 different 
6- and 12-atom pathways. That study required the use 
of an accurate and efficient interatomic potential model: 
in this tight-binding model reparameterized using 

all-electron density-functional calculations. 

After reparameterizing, we modify the functional form 
of tight-binding for small interatomic distances to over- 
come the collapse problem. This ensures that the poten- 
tial is suitable for phase transition studies and molecular 
dynamics simulations. The collapse problem for tight- 
binding models is caused by unphysically large overlap at 
small distances creating a low energy binding state; by 
modifying the functional form using short-range splin- 
ing, the collapse problem can be avoided. This paper 
provides the details of the model used in the previous 
phase transformation study of Ref. 3] and describes a 
general solution to the collapse problem. 

Tight-binding is a parameterized electronic structure 
method for calculation of total energies and atomic forces 
for arbitrary structures. It is an empirical model that 
can reproduce density-functional results for a range of 
structures yet requires orders of magnitude less compu- 
tational effort. The parameters of the model are deter- 
mined by fitting to a database and the range of appli- 
cability is determined by comparison to structures not 
in the database. The end result is a model that bal- 
ances three competing properties — efficiency, accuracy, 



and transferability — which make it applicable to a vari- 
ety of important structures. 

We fit our model to total energies and electron eigen- 
values for several crystal structures over a range of vol- 
umes to produce a transferable model for the study of the 
a^ui transformation^ Our fitting database is chosen to 
sample a large portion of the available phase space of 
parameters while constraining those parameters as much 
as possible. The resulting model accurately reproduces 
total energies, elastic constants, phonons, and point de- 
fects; all of which are necessary for transformation mod- 
eling. In addition, the functional forms are modified for 
small distances to overcome the unphysical collapse prob- 
lem; this is necessary for phase transitions and molecular 
dynamics which sample small interatomic distances. 

Section [D] describes tight-binding as a parameterized 
electronic structure method, the functional forms for ti- 
tanium, the modifications for short distances, our fitting 
database and our method of optimization. Section IIIII 
gives the optimized parameters, and tests our model 
against total energies, elastic constants, phonons, and 
point defect formation energies for a, u>, and bec Ti. The 
point defect formation energies are used to compare our 
parameters to those of Mehl and Papaconstantopoulos 4 
and Rudin et al.^ and to demonstrate the efficacy of our 
modification of the short-range Hamiltonian and overlap 
functions. 



II. METHODOLOGY 

A. Tight-binding formulation 

Electronic structure methods separate the total en- 
ergy of a crystal into an ionic contribution and an elec- 
tronic contribution derived as the solution to a Hamil- 
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tonian problem. Treating electrons as non-interacting 
fermionic quasiparticlcs permits an appropriate one- 
particle solution, 6 To numerically solve the electronic 
problem requires a set of basis functions <j>i, in terms 
of which the matrix Hij of the Hamiltonian operator and 
overlap matrix Sij are 

= {</>i\H\^j),Sij = (4>i\(t>j}- 
These matrices give the eigenvalue equation, 

Hlp n = tnSlpn, (1) 

where the electronic contribution to the total energy in- 
cludes the term 

2 J2 € " > 

e„<E F 

with Fermi energy Ef. The Hamiltonian contains in- 
formation about the wavefunction solutions themselves 
(e.g., density-functional theory). Typically, the wave- 
functions must be found self-consistently, which increases 
the computational requirements. 

In the tight-binding method, approximate Hamiltonian 
and overlap matrices are constructed by assuming atom- 
centered orbitals in a two-center approximation. This 
technique is related to the linear combination of atomic- 
like orbitals (LCAO) method, which uses a basis 0, of 
solutions to the isolated atomic Schrodinger equation up 
to some energy and angular momentum quantum num- 
bers (nl): (j) n i m {r) = f n i{\r\)Y lm (r/\r\). Tight-binding 
Hamiltonian and overlap functions are calculated inde- 
pendently of the local environment which increases effi- 
ciency but at the expense of transferability. 

Empirical tight-binding eliminates explicit basis func- 
tions from the problem and parameterizes the Hamilto- 
nian and overlap matrices in terms of simple two-center 
integrals^ The basis is chosen to be angular momentum 
solutions Im up to some maximum I value: For a maxi- 
mum I = 1 we use s, p x , p y and p z as the basis functions; 
for a maximum of I = 2, we add in the five rf-orbitals 
d X y, d yz , d zx , d x 2_ y 2 and d 3z 2_ r 2. The Hamiltonian and 
overlap matrices are written as sums of parameterized 
functions hi m ,l'm'(r) and s/ m ,;' m '(r) where r = R t - Rj 
is the separation between two atoms i and j. The two- 
center approximation allows these functions to be simpli- 
fied further according to the angular momentum compo- 
nents of the basisi For example, h PziPz (r) separates into 
two symmetrized integrals 

hp z , P z W = h ppa (r) cos 2 9 Z + h pplT {r) sin 2 Z , 

where Z is the angle between r and the z axis. The 
higher rotational angular momentum integral h pp s(r) is 
zero because a p orbital has a maximal azimuthal quan- 
tum number of 1 along the z axis. The integrals h ppa (r) 
and h ppn (r) are functions of only the distance of separa- 
tion r = \Rj — Rj'\. We write each Hamiltonian and 
overlap integral in these symmetrized functions; for a 



model with an spd basis, there are ten integrals (for h and 
s) to be determined: (ssa), {spa), (ppcr), (ppn), (sda), 
{pda), (pdir), (dda), (ddw), and (ddS). The Hamiltonian 
and overlap matrices are then computed for an arbitrary 
atomic arrangement. The total energy of the system is 
given by the eigenvalues of the Eq. £n and the ionic 
contribution: 

^total = 2 2^ £ " + ^(-Rnucloi), 

where V does not depend on the electronic states of the 
system. 

We use functional forms developed at the U.S. Naval 
Research Laboratory (NRL) that do not use an explicit 
external pair potential but instead has environment de- 
pendent onsite energies jSuLili The onsite Hamiltonian ele- 
ments e s , e p , and ed are not constants, but rather, depend 
on the distances of neighboring atoms to approximate 
three-body terms^ The onsite energies eij are functions 
of the "local density" pi with four parameters 

£M = a i + b lPi /3 + c iPi /3 + dlPi, (2) 

where 

Pi = cxp(-A% )/ c (r. y - ). (3) 
The smooth cutoff function / c (r) is 

'■M-( 1+ -»(Tr))"- < 4 > 

The intersite functions hu> m (r) and su> m (r) are given by 
three parameters each 

hwm(r) = {ew m + fli'mr) exp(-g? llm r)f c (r), 
swm{r) = (e Wm + fii'mr) exp(-g u , m r)f c (r). 

The squared parameters gwm and gw m guarantee the 
exponential terms to decay with increased distance. 

The overlap and Hamiltonian functions have an un- 
fortunate behavior for small distances r which can lead 
to catastrophic failure in the Hamiltonian problem. The 
functional form in Eq. (jSJ is exponentially damped as 
r grows; in reverse, this means that our intersite func- 
tions grow exponentially as r becomes small. As h or s 
between two atoms grow in magnitude they increase the 
bonding between the two respective atoms; as |s| — > 1 the 
energy of the bond grows as 1/(1 — |s|). When the bond 
energy grows, the bonding state is populated while the 
antibonding state is not; this results in a net attractive 
force between the two atoms. As the interatomic distance 
shrinks, the entire overlap matrix S ceases to be positive- 
definite, and the Hamiltonian problem of Eqn. Q is no 
longer solvable. This causes the "collapse problem" in 
molecular dynamics: two atoms come close to each other 
and see a large attractive force that pulls them towards 
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FIG. 1: Interpolated intersite function with short-range 
spline. The parameterized function h(r) grows exponentially 
as r approaches zero, though the function is only sampled in 
the fitting database down to flmin- At r = i? m in, we replace 
the function with a quartic spline that matches the value, first 
and second derivatives at -R m i n ; the dashed curve shows the 
growth of the original function. The spline smoothly goes to 
a constant value in a width of a. Only one adjustable pa- 
rameter, a, is added to the entire fitting database, as a is the 
same for all functions. 
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FIG. 2: Energy of Ti dimer calculated with tight-binding us- 
ing short-range splining. Without any short-range splining, 
the overlap matrix becomes artificially large, creating a bond- 
ing state with very low energy at small distances; at 1.92A, 
the tight-binding dimer Hamiltonian problem becomes un- 
solvable. As described in the text, by short-range splining of 
the Hamiltonian and overlap functions, the model is stable 
and becomes repulsive at small distances. 



each other until S is not positive-definite. In actuality, 
the Hamiltonian problem is not meaningful even before 
S is not positive-definite, because the model predicts a 
bond with an unphysically low energy. In a real mate- 
rial, the growth in bonding is counteracted by Coulumb 
repulsion: a two-electron term that is not included in the 
tight-binding formalism. 

Short-range splining. To resolve this, we modify the 
intersite functions to keep the overlap matrix S positive 
definite. Because our fitting database includes only in- 
teratomic distances larger than some minimum distance 
-Rmin, the functional form is guaranteed to be correct only 
for r > -Rmin- Below -R m in, we smoothly interpolate both 
hwm{r) and sw m (r) to a constant value. The interpola- 
tion is performed with a quartic spline, from r = -R m in 
down to r = -R m m — c; below i£ m i n — a, the function takes 
on a constant value. We choose spline values to enforce 
continuity of value and the first and second derivatives; 
the final functions for both hw m {r) and su> m (r) are 



{h(r) : r > i? min , 

frspiine(O) : r < i? min - a, 

h sp iinc(r - -Rmin + o) ■ otherwise, 

(6) 



where 

1 1 / 1 \ u 3 

h sp y mc (u) = h --o-h' + —o- 2 hg+ lah' - -a ftp J 

for u in [0, a], and ho, h' , and h' ' are the value, first, and 
second derivative of h(r) at i? m in- Figure^shows this in- 
terpolation schematically. While we smoothly interpolate 
har m and su> m , we retain the environment-dependent on- 
site terms; this has the effect of reducing the strength 
of bonding while the onsite energy continues to grow — 
effectively producing a pair repulsion between atoms at 
small distances. 

Figure [21 illustrates the collapse problem for the Ti 
dimer and how short-range splining stabilizes the model 
for small distances. As the distance between the two 
atoms decreases, a bonding state with an artificially low 
energy decreases the dimer energy. The precipitous drop 
in the energy of this bonding state is due to an increase in 
the overlap; at 1.92A, the overlap matrix becomes non- 
positive definite, and the eigenproblem is no longer solv- 
able. A a value of 0.5A or 0.25A makes the dimer stable; 
this is necessary but not sufficient to solve the collapse 
problem for all cases. 

Our parameterization has 74 parameters to be opti- 
mized, plus 3 fixed parameters. The cutoff function / c (r) 
has two fixed parameters Rq and Zq , while the minimum 
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distance R mm is set by the database. There are 10 Hamil- 
tonian and 10 overlap functions, each with 3 parameters 
for a total of 60 parameters. The 3 onsite energy func- 
tions have 4 parameters each, and a single parameter A 
for the density gives 13 parameters. Finally, the short- 
range spline range parameter a is determined using the 
dimer, and testing with molecular dynamic calculations 
and defect relaxations^ 



B. Fitting database 

We compile a database of electronic structure calcula- 
tions of several crystal structures using full-potential lin- 
earized augmented plane wave (FLAPW) calculations 11 
with the wien97 program suited We use the general- 
ized gradient approximation (GGA) for the exchange- 
correlation energy^ The sphere radius is Rmt = 
2.0 bohr = 1.06 A; there is a negligible charge leakage 
of 10~ 8 electrons. The planewave cutoff K max is given 
by i?MT^max = 9; this corresponds to an energy cut- 
off of 275 eV. The energy cutoff is not as large as re- 
quired in a typical pseudopotential calculation because 
the planewaves are only used in the interstitial regions 
away from atom centers. The charge density is expanded 
in a Fourier series; the largest magnitude vector in the ex- 
pansion G max is 18 bohr -1 (34 A -1 ). Local orbitals are 
used for the s, p, and d solutions inside the spheres^i 
Our core configuration is Mg with semi-core 3p states 
represented by the local p orbitals; our 4s, 3c?, and 4p 
states are the valence orbitals. A Fermi-Dirac smear- 
ing of 20 mRyd (272 meV) is used to calculate the total 
energy^* 

Table [Q shows a summary of the fitting database; it 
consists of the total energies and eigenvalues on a &-point 
grid for several crystal structures. Five structures are 
used: simple cubic (sc), body-centered cubic (bcc), face- 
centered cubic (fee), hexagonal closed-packed (a), and 
omega (to). The three cubic structures are calculated over 
a range of volumes, while the hexagonal structures are 
calculated only at FLAPW equilibrium volumes and c/a 
ratios. The eigenvalues in each structure are each shifted 
by a constant amount so that the sum of the occupied 
bands (smeared by 272 meV) is the total energy. We 
calculate the 9 lowest bands per atom above the semicore 
3p states; these represent the 4s, 3d and 4p states both 
below and above the Fermi level. We use the lowest 6 
bands at each /c-point for fitting the cubic structures, 9 
bands for a, and 12 bands for u>. 

In addition to eigenvalues on a regular grid, we include 
eigenvalues at high symmetry points and directions in the 
Brillouin zone to aid in fittingiiSiii For the three cubic 
structures, we calculate the eigenvalues at several high- 
symmetry points and directions (10 for bcc and sc, and 12 
for fee) and then decompose the electronic wavefunctions 
in terms of the symmetry character of the eigenvalues^ 18 
Again, we use the lowest 6 states for the high-symmetry 
points. We are careful not to fit too many eigenvalues 



TABLE I: Crystal structures used in tight-binding fitting 
database. Five different crystals structures are used, with 
five volumes for each of the cubic crystal structures. The 
lattice constant ao, volume per atom, and nearest neighbor 
distance and multiplicity for each structure is listed. The 
equilibrium lattice constant for each structure is marked with 
an asterisk. The same fc-point mesh is used for all volumes 
of a given structure, and is constructed using the prescription 
of Ref. [l4lT5j| . The smallest distance to appear in this fitting 
database is f? m i n = 2.350 A. 



Structure 


a (A) V/atom (A 3 ) 


nn (A) fc-point mesh 


bcc 


2.887 


12.03 


2.500 x8 shifted 5x5x5 




3.060 


14.32 


2.650 x8 (44 points) 




3.281* 


17.66 


2.841 x8 




3.406 


19.76 


2.950 x8 




3.579 


22.93 


3.100 x8 


fee 


3.747 


13.16 


2.650 x 12 unshifted 5x5x5 




3.960 


15.52 


2.800 x 12 (47 points) 




4.127* 


17.57 


2.919 xl2 




4.384 


21.06 


3.100 xl2 




4.596 


24.27 


3.250 xl2 


sc 


2.350 


12.98 


2.350 x6 shifted 5x5x5 




2.500 


15.62 


2.500 x6 (35 points) 




2.645* 


18.50 


2.645 x6 




2.800 


21.95 


2.800 x6 




2.950 


25.67 


2.950 x6 


a 


2.952* 


17.69 


2.952 x 12 unshifted 5x5x2 




(c/a = 


1.588) 


(42 points) 




4.600* 


17.23 


2.656 x3 unshifted 3x3x4 




(c/a = 


0.613) 


(35 points) 



FLAPW equilibrium lattice constant 



at high-symmetry points, since the lowest 9 bands in the 
GGA band structure may not correspond to those pre- 
dicted in our spd basis^ 

Because our fit includes the electron eigenvalues, we 
expect our model to reproduce both total energies and 
energy derivatives. Phonons and elastic constants can be 
written in terms of the forces on atoms due to small dis- 
placements; the Hellman-Feynman theorem relates the 
force on an atom Ri to the eigenvalues as 

e n <E F 1 

Thus, the electron eigenvalues of the bulk crystal contain 
information about phonons and elastic constants. 

C. Optimization of parameters 

The parameters are optimized to minimize the mean 
squared error. We use the non-linear least-squares mini- 
mization method of Levenberg-Marquardt with a numer- 
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ical JacobianiiS We weight each /c-point by unity, and the 
resulting total energy by 200; accordingly the total ener- 
gies are weighted approximately the same as the fc-point 
data. We initialize our parameters usin g th e Hamilto- 
nian and overlap values for Ti from Ref. 16] adapted to 
our functional form. We then fit only the environment- 
dependent onsite terms to the band-structure of the cubic 
elements. After an initial fit is found, we include the hop- 
ping terms in the optimization. We proceed using only 
the cubic band structure, then the cubic band structure 
and total energies, and finally all structures and energies. 
After a new minimum is found, we check each function 
to see if the minimization has made the exponential term 
gtll'm) t 00 large; this corresponds to making the entire 
function approximately zero over the sampled range of 
r values. We remedy this by resetting the e, /, and g 
parameters to 0, 0, and 0.5. Several fitting runs are per- 
formed until the entire fit set is accurately reproduced. 

III. RESULTS 
A. Parameters and fitting residuals 

Table [H] lists the parameters of the optimized tight- 
binding model. Figure shows the hopping integrals 
h(n> m )(r) and s (H / ro )(r) for a range of volumes; the R min 
in the database is 2.35 A, and we interpolate each func- 
tion to a constant value below i? m i n . Finally, Figure 0] 
shows the environment-dependent onsite energies as a 
function of volume for an hep crystal with c/a — 1.588. 

To use the potential for phase-transformation studies, 
a was determined by testing the stability of (1) the dimcr, 
(2) molecular dynamics runs, and (3) defect relaxations. 
While the lowest energy pathways studied by Ref. Q 
have distances of closest approach of 2.6 A, there were 
possible pathways where atoms approached within 2.3 A 
of each other. Without short-range splining, calculations 
of energies of structures with distances below our i? m ; n 
value can become problematic. Initially, a a value of 
0.529 A was chosen based on the dimer; however, defect 
relaxation calculations showed that a value of 0.265 A 
was necessary to ensure stability for some of the point 
defects. 

Table IIIII lists the errors in our tight-binding model 
with respect to the fitting database. Our average to- 
tal energy errors are approximately 1 meV; root-mean 
square errors in the fc-point energies are approximately 
100 meV. The tight-binding parameterization adequately 
reproduces the database energetics. To test transferabil- 
ity, we compare to properties outside of this database. 

B. Total energies 

FigurcElshows the tight-binding total energy as a func- 
tion of volume for a and u>. These curves were not 
included in the fitting database; only the two points 



TABLE II: Tight-binding parameterization for titanium. The 
onsite parameters are given for the s, p, and d orbitals. Each 
term is density dependent; the parameter in the density de- 
pendence is A. The cutoff function has fixed parameters Ro 
and lo- Next, the intersite Hamiltonian and overlap elements 
are given for each of the 10 symmetrized (11' m) combinations. 
Below -R m j n , each intersite function is smoothly interpolated 
to a constant value over the range a. 

a, (cV) b, (cV) Cl (cV) d, (cV) 

s: -3.272 x 10° 3.714 x 10 2 8.029 x 10 3 7.879 x 10 4 



p: 


4.974 x 10° 3.747 x 


10 1 -1.874 x 10 3 2.721 x 10 4 


d: 


3.632 x 10 _1 3.238 x 


10 1 8.877 x 10 1 9.355 x 10 2 




— a L + bip t ' + 


Cip/ +dip i 







Pi — 


cxp{~X 2 r ij )f c (r ij ) : A 2 = (0.3620 A) 1 J3 


f c (r) 


= ( 


1+cxp (Sr) 


\ Ho — 6.615 - 
/ to — 0.2b4o 


K ® 






e tU'ml 


/ril'ml 




sscr: 


h— 


-1.086 x 10 2 cV 


-3.900 x 10 3 oV/A 


0.3277 A 




s— 


9.277 


-2.624 A -1 


0.8357 A 


spa: 


h= 


-1.793 x 10 3 eV 


8.066 x 10 2 cV/A 


0.4926 A 




s— 


-11.81 


0.02523 A" 1 


0.5993 A 


ppa: 


h= 


-4.865 x 10 2 cV 


1.816 x 10 2 cV/A 


0.6929 A 




s— 


0.08093 


-1.351 A -1 


1.036 A 


ppn: 


h= 


1.202 X 10 1 cV 


-8.252 x 10° cV/A 


0.8925 A 




s— 


4.478 


-0.2899 A -1 


0.8026 A 


sda: 


h= 


-5.537 x 10 2 cV 


3.096 x 10 2 oV/A 


0.4772 A 




s— 


-4.331 


-5.085 A -1 


0.4498 A 


pd<r: 


h= 


-2.338 x 10 2 oV 


9.994 x 10 1 oV/A 


0.6321 A 




s— 


0.02557 


-3.383 A -1 


0.5728 A 


pdrr: 


h= 


-4.979 x 10° cV 


7.855 x 10 _1 cV/A 


1.617 A 




s— 


0.1943 


2.308 A -1 


0.5882 A 


dda: 


h= 


1.706 X 10 2 eV 


-1.150 x 10 2 oV/A 


0.5266 A 




s= 


-0.9905 


0.7605 A -1 


0.7990 A 


ddrr: 


h= 


9.920 X 10° oV 


3.538 x 10 1 cV/A 


0.5366 A 




s— 


-1.490 


-1.498 A" 1 


0.5213 A 


ddS: 


h= 


1.109 x 10 3 cV 


-6.205 x 10 2 cV/A 


0.3340 A 




s— 


15.58 


-5.276 A" 1 


0.4412 A 



0> s }(!I'm)( r ) = ( e (!!' m ) + f(ll m ') r ) cx P(-9^;' m ) r )/<=( r ) El 

fl min = 2.350 A, g = 0.265 A ® 



indicated. We reproduce both the slightly lower en- 
ergy of u) over a predicted by pseudopotential methods 3 
and FLAPW calculations, as well as the slightly lower 
equilibrium volume of to. The three cubic structures 
were included in the fit and have errors on the order of 
3 meV/atom (c.f., Table 1111}) . This shows a wide range 
of applicability for our model under pressure. 

C. Elastic constants and phonons 

Table IIVI shows the equilibrium lattice constants and 
elastic constants for a, lu, and bec for our tight-binding 
model. The GGA numbers correspond to the elastic con- 
stants found using VASPi^L^i Elastic constant combina- 
tions which do not break symmetry such as C\\ + C12, 
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FIG. 3: Tight-binding intersite Hamiltonian and overlap func- 
tions. The parameterized hopping integrals are shown for dis- 
tances from 2.4 A to 3.5 A. The -Rmm in the fit is 2.35 A; below 
this, these functions are smoothly interpolated to a constant 
value. Circles represent a integrals, squares n integrals, and 
triangles 5 integrals; black is for s, dark gray for p, and light 
gray for d. 



TABLE III: Fitting errors in total energy and fc-points for 
tight-binding model. For each structure, we report the abso- 
lute error in the total energy (first line) and the RMS error 
in all fc-point energies in the fit set (second line). The to- 
tal energy errors are on the order of 1 meV, while the RMS 
band-structure errors are on the order of 100 meV. 





low volume 


equilibrium 


high volume 


bcc 


1.64 meV 


0.957 meV 


4.31 meV 




200 meV 


104 meV 


110 meV 


fee 


-1.79 meV 


1.25 meV 


-0.821 meV 




136 meV 


87.1 meV 


114 meV 


sc 


-0.0190 meV 


-0.115 meV 


-1.60 meV 




435 meV 


195 meV 


140 meV 


a 




-1.66 meV 
69.1 meV 








-0.00993 meV 
67.9 meV 






volume/atom (A ) 



FIG. 4: Tight-binding onsite energy terms for hep structure. 
The onsite energies are environment dependent in our model; 
we show the variation with respect to the volume of an hep 
crystal with c/a — 1.588. The low volume of 10 A 3 has a 
lattice constant of 2.44 A, and the high volume of 25 A 3 has 
a lattice constant of 3.31 A. The equilibrium hep volume is 
17.56 A 3 . 
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FIG. 5: Comparison of tight-binding energy as a function of 
volume for a and u with first principles data. The two filled 
points were included in the fit; the lines are FLAPW total 
energies. Our tight-binding model reproduces the fit data — 
slightly lower ground-state energy and equilibrium volume for 
uj — and the equation of state of the full-potential calculations. 



TABLE IV: Lattice parameters and elastic constants in units 
of GPa for a, ui, and bcc Ti from tight-binding, GGA, and 
experiment. GGA corresponds to the elastic constants found 
using vasp^"*— The experimental a elastic constants are mea- 
sured at 4K,—- and the bcc elastic constants at 1238K, ?3 
Our tight-binding model reproduces the GGA elastic con- 
stant combinations that preserve the symmetry of the struc- 
ture (e.g., Cii + C12), but has larger error with those that 
break it (e.g., C44). The deviation between the bcc exper- 
imental elastic constants and our calculations is due to the 
high temperature needed to stabilize the bcc structure in Ti. 



a (A) c (A) Cn gig C13 C33 C44 
Tight-binding 



a 


2.94 


4.71 


155 


91 


79 


173 


65 




4.58 


2.84 


184 


90 


52 


261 


100 


bcc 


3.27 




87 


112 

GGA 






31 


a 


2.95 


4.68 


172 


82 


75 


190 


45 




4.59 


2.84 


194 


81 


54 


245 


54 


bcc 


3.26 




95 


110 






42 



Experiment 

a 2.95 4.68 176 87 68 191 51 

bcc 3.31 134 110 36 



C13, C33 in the hexagonal crystals, and Cn +2C12 in bcc 
are reproduced within approximately 10%. However, the 
symmetry breaking elastic constant combinations such as 
Cn — C12 and C44 have larger errors. It is worth not- 
ing that none of this data, except for the bulk modulus 
of bcc, appears in any form in the fitting database; the 
agreement is a consequence of reproducing the electron 
eigenvalues. 

We calculate phonons using the direct-force 
methodi22i2£i2§i2I We calculate the forces on all atoms 
in a supercell where one atom at the origin is displaced 
by a small amount. The numerical derivative of the 
forces with respect to the displacement distance approx- 
imates the force constants folded with the translational 
symmetry of the supercells. The Fourier transform of 
the force constants gives the dynamical matrix, and 
its eigenvalues give the phonon frequencies^ For q 
vectors commensurate with the supercell, the phonon 
frequencies are exact; for incommensurate q vectors, the 
calculated phonons are a Fourier interpolation between 
exact values. Our supercells are 4x4x3 for a, 3x3x4 
for u>, and 4x4x4 simple cubic cell for bcc; in all cases, 
a 2 x 2 x 2 fc-point mesh is used in the supercell. 

Figures H3 and [S] are the predicted phonon dis- 
persions for our tight-binding model, calculated at the 
equilibrium volumes for each structure. The a phonons 
match the experimental values well for the high energy 
phonons optical and acoustic branches; these are impor- 
tant for modeling the shuffle during martensitic trans- 
formation. The deviation from experiment for small q 
corresponds to our mismatch in elastic constants. The u> 



7 



Koo] [%o] [ooy 




FIG. 6: Comparison of tight-binding phonons for the a phase 
with experimental phonon data. The crosses are the experi- 
mental phonon frequencies at 295Ki— ■ The deviation from the 
experimental values at small q corresponds to the mismatch 
in the a elastic constants. Our tight-binding model does well 
for the high-energy optical and acoustic branches which are 
important for modeling the a— >lu transformation. 
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FIG. 7: Predicted u) phonons from tight-binding. As expected 
from the c/a ratio of 0.620, the phonon modes are stiffer along 
the [00£] direction than the basal plane directions [£00] and 
[0*0]. 
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FIG. 8: Predicted bcc phonons from tight-binding. At T = 0, 
the bcc phase in Ti is unstable, as shown by the imaginary 
phonon frequencies. The dip in the branch is near the L- 
§[111] phonon, which corresponds to the bcc— >oj transforma- 
tion pathway. The imaginary phonon for T-[011] corresponds 
to the bcc— >a transformation mechanismi— 

phonons are expectedly stiffer along the c axis than in the 
basal plane due to the low c/a ratio. The bcc phonons 
show phonon instabilities corresponding to the bcc— >lu 
transformation (L-|[lll] phonon) and the bcc— >a trans- 
formation (T-[011] branch). 30 

D. Point defects 

Table IVl shows the formation energies of point defects 
for a and u> at the equilibrium volumes for our tight- 
binding model. All a calculations are performed with a 
4x4x3 (96 atom) supercell and all u with a 3 x 3 x 4 (108 
atom) supercell. No point defect information is included 
in the initial fit; we reproduce the GGA formation ener- 
gies for all of the point defects considered. This indicates 
that our tight-binding model is applicable to the study 
of the a— >w transformation path, where atoms move out 
of their equilibrium configurations and often close to one 
another. 

The formation energies of point defects shows some im- 
provement of our model over two existing modelsA^ The 
potential by Rudin et al. uses the same functional forms 
as our potential without short-range splining for hopping 
and overlap functions; Mehl and Papaconstantopoulos 
use the same onsite function form, but adds additional 
quadratic parameters to the hopping and overlap func- 
tions in Eqn. ©. All three potentials use the same onsite 
functional forms. For all three potentials, the binding en- 



S 



TABLE V: Point defect energies in eV for a and to Ti from 
tight-binding for different parameterizations. GGA refers to 
the defect formation energies calculated with VASP; TB to 
the parameterization in this work; NRL to the parameteriza- 
tion by Mehl and Papaconstantopoulos)- and LANL to the 
parameterization by Rudin et al. 5 NN refers to the distance 
of closest approach for two atoms in each defect in TB. The 
formation energies are calculated after relaxation. The de- 
fects marked coll. fell victim to the "collapse problem" dur- 
ing relaxation. The a-tetrahedral site is unstable, relaxing to 
form a dumbbell along the [0001] direction in GGA. The u>- 
hexahedral site is very close to the w-tetrahedral site*— Many 
of the interstitial defects sample small distances, requiring the 
use of short-range splining to stabilize the defects. 



Defect 


GGA 

n 


TB 

defects 


NRL 


LANL 


NN [Aj 


Octahedral 


2.58 


2.89 


1.31 


2.55 


2.50 


A 


Tetrahedral 






unstable 






Dumbbell-[0001] 


2.87 


2.81 


1.81 


coll. 


2.18 


A 


Vacancy 


2.03 


1.88 


1.51 


1.92 


2.83 


A 


Divacancy-AB 


3.92 


3.83 


3.73 


3.68 


2.81 


A 




i,i 


defects 










Octahedral 


3.76 


4.11 


3.20 


3.67 


2.30 


A 


Tetrahedral 


3.50 


3.58 


2.86 


coll. 


2.21 


A 


Hexahedral 


3.49 


3.86 


2.88 


4.37 


2.28 


A 


Vacancy-A 


2.92 


2.85 


2.99 


3.25 


2.60 


A 


Vacancy-B 


1.57 


1.34 


1.01 


1.90 


2.62 


A 



ergies versus volume, elastic constants, and phonons are 
similar, though Rudin's more accurately captures the low 
frequency a phonons. However, point defect formation 
energies are better predicted by our tight-binding param- 
eterization. 



The short distances sampled by the point defects em- 
phasize the need for short-range splining of both the over- 
lap and Hamiltonian functions. The collapse of two de- 
fects in Rudin et a/.'s model is due to the growth of 
the overlap matrices; the lower energies predicted by 
Mehl and Papaconstantopoulus could be due to overly 
large overlap elements at short-distances as well. Inter- 
stitial defects, like phase transformation pathways, can 
sample interatomic distances smaller than the smallest 
distance included in the fitting database; without short- 
range splining, this can lead to artificially lower energies, 
or even collapse. Without short-range splining, all three 
tight-binding parameterizations fail for the Ti dimer at 
small distances: 1.92 A for this work, 1.76A for Mehl and 
Papaconstantopoulos, and 1.28 A for Rudin et al. The use 
of short-range splines provides a solution to the collapse 
problem for non-orthogonal tight-binding models. 
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IV. CONCLUSION 

We present an accurate and transferable tight-binding 
model with parameters determined by density-functional 
calculations. It reproduces structural energies with pres- 
sure, elastic constants, phonons, and point defect ener- 
gies. By fixing the short-range behavior of the potential, 
point defects can be accurately computed, which allows 
the calculation of energy barriers for phase transforma- 
tion pathways. The wide range of applicability makes it 
particularly well suited to the study of martensitic phase 
transformations, such as a— and short-range splines 



represent a solution to the potential collapse problem of 
non-orthogonal tight-binding models. 
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